G x E expression

The question is how much G x E we see in the RNAseq data.

To get started on this I ran the file “ProcessCount.R”. For this the data was TMM normalized, voom transformed, and then a gtXenvironment model was fit in limma.

library(limma)
library(tidyverse)

── Attaching packages ────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 2.2.1     ✔ purrr   0.2.4
✔ tibble  1.3.4     ✔ dplyr   0.7.4
✔ tidyr   0.7.2     ✔ stringr 1.2.0
✔ readr   1.1.1     ✔ forcats 0.2.0
── Conflicts ───────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(stringr)
library(modelr)
library(broom)

Attaching package: ‘broom’

The following object is masked from ‘package:modelr’:

    bootstrap

library(multidplyr)
library(VCA)

Get the limma objects

load("voom.int.Rdata")

How many genes show GxE?

generate vector of interaction coefficients and then topTable

colnames(fit.int$coefficients)
int.coef.n <- colnames(fit.int$coefficients) %>% str_which(":") 
topTable(fit.int,coef = int.coef.n)

Not working

let’s just do an ANOVA for each gene

counts <- as.tibble(counts.voom.gr$E) %>% 
  mutate(geneIndex=1:n())
counts[1:10,1:10]

counts.long <- counts %>% gather(key="name",value="expression", -geneIndex) %>%
  mutate(gt=factor(str_extract(name,"RIL_[0-9]*")), 
         env=factor(str_extract(name,"CR|UN")), 
         rep = factor(str_extract(name,"Rep[0-9]"))) %>%
  select(-name) %>%
  group_by(geneIndex) %>% nest()
counts.long

counts.long$data[[1]] 

#test.df <- counts.long$data[[1]] 

Now fit an anova and do a VCA for each gene in parallel

fitAnova <- function(df) aov(expression ~ gt*env,data=df)

getPvals <- function(fit) {
  pvals <- summary(fit)[[1]]$`Pr(>F)`
  names(pvals) <- dimnames(summary(fit)[[1]])[[1]] %>% str_trim()
  data.frame(t(pvals))
}

fitVCA <- function(df) {
  tmpvca <- anovaVCA(expression ~ gt*env,Data=as.data.frame(df))
  data.frame(t(tmpvca$aov.tab[,"%Total"]))
}

cl <- create_cluster(cores=4)
counts.long <- partition(counts.long,geneIndex,cluster = cl)
cluster_library(cl,c("purrr","broom","VCA","stringr"))
cluster_copy(cl,fitAnova)
cluster_copy(cl,getPvals)
cluster_copy(cl,fitVCA)

system.time(
  counts.results <- counts.long %>% mutate(aov = map(data,fitAnova), 
                                        glance = map(aov,glance),
                                        pvals = map(aov,getPvals),
                                        varcomp = map(data,fitVCA)) %>%
    select(-data, -aov)
)
system.time(cluster_rm(cl,"counts.long"))
rm(counts.long)

system.time(counts.results <- collect(counts.results))
save(counts.results,file="GxE_models.Rdata")

(aovsummary <- counts.results %>% unnest(glance,.drop = TRUE))

(pvals <- counts.results %>% unnest(pvals,.drop = TRUE))

(varcomp <- counts.results %>% unnest(varcomp,.drop = TRUE))

PVals

Number of genes (out of 24619) that have a FDR adjusted pvalue of less than 0.05:

pvals %>% ungroup() %>%
  summarise(gt=sum(p.adjust(gt,method="fdr") < 0.05),
            env=sum(p.adjust(env) < 0.05),
            gt.env=sum(p.adjust(gt.env,method="fdr") < 0.05))

Variance components

Histograms of % variance explained

pl <- varcomp %>%
  ungroup() %>%
  select(-total) %>%
  gather(key = "component", value = "percent.variance",-geneIndex) %>%
  ggplot(aes(x=percent.variance)) +
  facet_wrap(~ component, scale="free_y") +
  geom_histogram()
pl

focus in on gt x env histogram, non-zeros

pl <- varcomp %>%
  ungroup %>%
  filter(gt.env > 0) %>%
  ggplot(aes(x=gt.env)) +
  geom_histogram(binwidth = 5)
pl

Number of genes with gt.env explaining X % variance:

varcomp$gt.env %>%
  cut(breaks=seq(0,35,5),right = FALSE) %>% table

.
  [0,5)  [5,10) [10,15) [15,20) [20,25) [25,30) [30,35) 
  23311    1121     151      29       5       1       1